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We discuss the dynamics of phase separation through the process of spinodal decomposition 
in a Fermi superfiuid with population imbalance. We discuss this instability first in terms of a 
phenomenological Landau theory. Working within the mean-field description at zero temperature, 
we then find the spinodal region in the phase diagram of polarisation versus interaction strength, 
and the spectrum of unstable modes in this region. After a quench, the spinodal decomposition 
starts from the Sarma state, which is a minimum of the free energy with respect to the order 
parameter at fixed density and polarisation and a maximum at fixed chemical potentials. The 
possibility of observing non-trivial domain structures in current experiments with trapped atomic 
gases is discussed. 
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The ordering of matter into different phases is a central 
preoccupation of many areas of physics, from condensed 
matter to cosmology. Hand in hand with the existence 
of different phases goes the question of the dynamical 
processes responsible for their formation, which may be 
equally important in determining what is observed in a 
given situation. Recent experimental advances in the cre- 
ation of degenerate atomic gases have begun to realize the 
prospect of a rich variety of new phases in atomic mat- 
ter, involving the hyperfine degrees of freedom, mixtures 
of different species, or spatial order on optical lattices. 
With each new phase comes the dynamical issue of how 
that phase will appear under laboratory conditions. One 
advantage offered by atomic systems is that the charac- 
teristic timescale h/ksT at nanokelvin temperatures is 
in the convenient millisecond range. 

The possibility of tuning interparticle interactions in 
a controlled manner has proven to be of particular 
significance lately. Magnetically tuned Fcshbach reso- 
nances have permitted the experimental investigation of 
the crossover from a Bose-Einstein condensate (BEC) 
of diatomic molecules to the Bardeen-Cooper-Schrieffcr 
(BCS) limit of weakly-bound Cooper pairs of fermionic 
atomsi2^££. 

It appears that when the numbers of atoms of the 
two species undergoing pairing are equal, the system 
forms a condensate with smoothly varying properties 
at low temperatures. With unequal numbers (we will 
call such a system 'polarised') there is the possibility 
of phase separation into a superfiuid of low polarisation 
(favoured by pairing) and a normal fluid of higher polar- 
isationi&^ii^M. 

With the occurrence of phase separation in these sys- 
tems now established, it is crucial to examine the dy- 
namics of this process, and that is the purpose of this 
Letter. Working within a mean-field approximation, we 
find that the dome of phase separation at temperatures 
below the tricritical poin t 14 > 15 contains a spinodal region 
where phase separation proceeds via a linear instability. 
In many recent papers this region is incorrectly identified 
with the phase coexistence regio n 16 i 17 i 18 i 19 (see Fig. []}. 




FIG. 1: (Colour online) Zero temperature mean-field phase 
diagram for magnetisation m/n versus interaction strength 
1/kFa. The phase separated region (PS) can be decomposed 
into two regions — the spinodal or unstable region and the 
metastable region (darker shade) — divided by the spinodal 
normal (sp-N) line (dashed) and the spinodal superfiuid (sp- 
SF) line (dot-dashed). In addition, the PS spinodal phase is 
divided in half by the line where the superfiuid density Q is 
zero (thin solid). The tricritical point (orange circle) is at 
m/n — 1 and 1/kFa ~ 2.37. Inset: schematic finite tempera- 
ture phase diagram of T c /ef versus m/n for a fixed interaction 
strength 1/kpa. Arrows indicate possible quenches into the 
spinodal region. 



We discuss the dynamics initiated by a quench into the 
spinodal region by finding the unstable density modes 
associated with the endpoint of the quench. The modes 
with the fastest growth rate give the characteristic size 
of the resulting domains of superfiuid and normal fluid. 
Though our description of these processes will not be 
quantitatively correct in the crossover region where the 
system is strongly interacting, we emphasise that this 
behaviour is a generic feature of systems possessing this 
type of phase diagram. In particular, our system shares 
many similarities with the problem of 3 Hc- 4 He mixtures 
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(see e.gi 20 ' 21 ). We will confine ourselves to the early 
stages of spinodal decomposition characterized by the ex- 
ponential growth of unstable modes. The emergence of 
a coarsening regime at later times, where domains scale 
with the time since the quench 2 ^, is a fascinating possi- 
bility that we will leave for future work. 

The rest of this paper is organized as follows. In 
the next section we introduce a simple phcnomenologi- 
cal model for the polarized system. We use this model 
to show how, with a simplifying assumption, mean-field 
calculations within the grand canonical ensemble may be 
applied to the early stages of spinodal decomposition. By 
adding dynamical assumptions, the equation satisfied by 
the sound velocity can be inferred. In section [II] we pro- 
vide a microscopic mean-field calculation of the spectra 
of unstable modes, which are discussed in some detail in 
section UTTl before we conclude. 

I. PHENOMENOLOGICAL MODEL 

We begin by discussing the phase diagram in phe- 
nomenological terms starting from a model free energy 
depending on the supcrfluid order parameter A and the 
density difference m = n.f — tit, following a similar ap- 
proach to the 3 He- 4 He system 2 ^ 

/ m (A,m)^|A| 2 + M |A| 4 + «|A| 6 

+ \x- 1 m 2 + im \A\ 2 . (1) 

The potential for A is the simplest one that can describe 
a first order transition. We are interested in r < so 
that for Tji — there is always superfluid order A ^ 
0. Further, the coupling to m has the obvious physical 
meaning that larger polarisations discourage pairing. We 
ignore for the moment additional couplings to the total 
density n = n-j- + 7i^— . 

Minimising ([1]) on the order parameter, gives the po- 
tential 

F(m) = mm/ m (A,m) . (2) 

Where the transition is first order, the phase coexistence 
region (m cr .gp, to c] ._n) is obtained by the usual tangent 
construction 2 ^. Inside the coexistence region one identi- 
fies a spinodal region (m sp _N, w sp _sf), where the suscep- 
tibility d^ n F(m) < 0, and phase separation proceeds via 
the growth of unstable modes. In general this represents 
an extremely difficult dynamical problem. We will make 
the simplifying assumption that following a quench into 
the spinodal region, the order parameter relaxes rapidly 
do its minimum at fixed m, while m(r), being a con- 
served quantity, begins to develop inhomogeneities on a 
much longer timescale. For the case of trapped gases, we 
will review the validity of this approach a posteriori. 

Often it is convenient, particularly in many body cal- 
culations, to work instead with the grand canonical po- 



tential /ft (A, h), obtained from f m in the usual way 
/ft(A, h) = min [/ m (A, m) - hm] 

m 

= lr\A\*+u\A\ i + v\Af-± Xn h\ (3) 

where r = r + 2jhx n and u = u — \xnl 2 < 0. If we 
assume u > 0, then for small 7 a second order transi- 
tion occurs when r changes sign. Increasing 7 causes the 
transition to become first order as u becomes negative. 
Minimising /ft(A,/i) on A gives the thermodynamic free 
energy Q(h) = minA /ft (A, h). At some critical h = h CT 
there is a discontinuity in the derivative (see Fig. [2]). In a 
uniform phase m = — 9ftSl so the boundaries of the phase 
coexistence region are m cr _N,cr-SF = — 9ftOL±. The two 
phases correspond to the two minima of fh(A,h): one 
at A = 0, the normal phase, and the supcrfluid phase at 
finite A. One can continue past /i cr on the metastablc 
minimum, rather than the true minimum. Such states 
are linearly stable, with a positive susceptibility — <9 2 f2 
(= [d^ n F(m)]~ 1 ), since one may easily see that 

a 2 n_a 2 /ft (d*f h y x (dm\ 2 

dh 2 dh 2 \dA 2 J \dA J ' U 

and the first term on the right hand side is — Xn < — 
the normal state susceptibility is assumed positive. Since 
the curvature of fl(h) remains negative, m takes values 
inside the coexistence region but the system remains uni- 
form. This is of course not the equilibrium state of the 
system, but the other phase must nucleate in order for 
phase separation to occur. 

This situation changes when the metastable minima 
merge with the maximum, for h corresponding respec- 
tively to m sp _N and to sp _sf for the metastable super- 
fluid and normal phases. As already discussed, inside 
the spinodal region we seek the constrained minimum of 
/ m (A,m) with respect to A at fixed m, and the spec- 
trum of unstable modes about this point. How should 
this programme be implemented for a many-body calcu- 
lation that provides instead the grand canonical potential 
/ft (A, h)l Fortunately, one may easily show that the sta- 
tionary points of / m (A,m) with respect to A coincide 
with those of /ft (A, h) at corresponding values of h and 
771. This is because f m (A,m) may be written as 

/ m (A,m) - f h (A,h*(A,m)) + h*(A,m)m , 

where h*(A,m) is the value of h that maximises 
fh(A,h) + hm. Thus m = —dh£l\h=h» ■ We can then 
easily see that the conditions, dAf m (A,m) = and 
0A/ft(A, h) = are equivalent when h = h*(A, m). 

In the spinodal region the unstable constrained mini- 
mum corresponds to a maximum of /ft, as may be seen 
from Eq. Q: d 2 n/dh 2 > requires d 2 f h /dA 2 < 0. 
In the context of the mean-field theory for the paired 
fcrmion system, this corresponds to the solution of the 
self-consistent equation in a magnetic field discovered by 
Sarma2£. 
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FIG. 2: (Colour online) Mean-field free energy Q(h/[i) for the 
Hamiltonian Q versus h/fj, at unitarity 1/kpa = (fi > 0) 
obtained evaluating the free energy density /h(A//z, /i//x) re- 
spectively at the local minimum for A/fi ~ 1.15 (green), at 
A = (blue) and at the local maxima (red), i.e. the Sarma 
state. The value of hj \i corresponding to the spinodal-normal 
point (sp-N), the spinodal-superfluid point (sp-SF) and the 
critical value for phase separation (cr) are indicated. At uni- 
tarity, the corresponding values of the magnetisation are re- 
spectively m cr _sF/n = JTisp-SF/n = 0, m sp .^/n ~ 0.73 and 
m cr _N/n — 0.93 (see Fig. [1]). Inset: Plots of the free en- 
ergy density /h(A//i, h/fi) for the values of h corresponding 
to /i S p-N, fosp-SF and h CI . Note that the thermodynamic free 
energy corresponding to the Sarma state has a positive cur- 
vature indicating instability. The cusp structure shrinks to a 
point at the tricritical point. 



m as specified in the model Eq. (JTJ) . m is written in terms 
of the distribution function of the majority quasiparticle 
distribution function (at zero temperature there are no 
minority quasiparticles) 



m(r) 



p 



™t(p,r), 



which obeys the Boltzmann equation 

[d t + v F • V r - 2 7 A V r <SA L • V p ] n T (p, r, t) = 0. 

where vp = vpp is the Fermi velocity. The linearized 
solution is expressed as 



"t(p,q,w) 



q • v F 



-2 7 A (5A L (q,w)<5(ep -fx) 



w — vp ■ q 
m(q,u>) = 2jA v(n)L(u)/\q\)5A L (c[,Lj) 

where L(x) = § log — 1 is the Lindhard function, and 
v{fx) the Fermi surface density of states. The dispersion 
relation of the linearized modes is then given by a solution 
of 



A iBio 2 7 A 

-Ww Qq 2 - Ru 2 

2 7 Aoi/(/i)i(w/|q|) -1 



= 0, 



(A + 4j 2 A 2 uL(c s )) (Q - Rc 2 ) - B 2 c 2 s = 0. (6) 

We will see that an equation of the same form emerges 
from the microscopic analysis of the next section, where 
the solutions will be further analyzed. 



These results are readily adapted to the presence of a 
trap potential V(r) using the following simple approxi- 
mation: At each point in space we find the constrained 
minimum of the order parameter for fixed local magneti- 
sation m(r)/n(r) corresponding to the density profiles 
before the quench. Spinodal decomposition therefore oc- 
curs where the local magnetisation lies in the spinodal 
region of the homogeneous phase diagram. Moving out 
from the centre of the trap corresponds to a vertical tra- 
jectory in Fig. [T] which will be displaced horizontally by 
a quench in 1/kpa. 

To make contact with the microscopic analysis of the 
next section, let us add appropriate dynamics to the 
above model. Expanding Ao + SA around some value we 
write a phenomcnological quadratic Lagrangian for the 
longitudinal (amplitude) and transverse (phase) modes 

C = B5A L SA T + |(<5A T ) 2 - I (V<5A T ) 2 



A 



(6A L ) 2 - 2 7 mA <5A L . (5) 



Eq. ([5]) includes all terms up to second order in the deriva- 
tives except for a (5 At) 2 term that will not change the 
sound velocity. 5Ap is coupled to the density difference 



II. MICROSCOPIC CALCULATION 

We turn now to the analysis of the microscopic problem 
described by the Hamiltonian 



H— fx„n 

: 9 



- ^ ^ ( e k Ma) c k<7 C : 



E4- 



y -k+q/2T l '-k+q/21 L --k'+q/21 C k'+q/2 
k,k',q 

where = k 2 /2m (we set H = 1) and where the scatter- 
ing length is introduced in the usual way: 



T . (7) 



lml x ^ 1 
g Ana V ^-j 2e k 



(8) 



The condition d 2 fh/dA 2 = 0, corresponding to a diver- 
gent susceptibility, was used to obtain the spinodal lines 
in Fig.[TJ As explained before, the unstable modes are to 
be found from the matrix response function (dynamical 
susceptibility) : 

X(r-r',t-t') 

/([n(r,t),n(r',0]> ([n(r, t), m(r>, t')])\ 
1 ([m(r,t),n(r',t')]) ([n(r, t), m(r', t')]) J ' W 
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Finding the spectrum of collective modes requires us to 
solve the equation 



detx- 1 (q,e(q))=0, 



(10) 



which defines the dispersion relation e(q). The unstable 
modes correspond to 3e(q) > 0. In practice, the response 
matrix ^ can be found making use of a path integral 
formulation, by expanding the action 



S[A lfl ,h] 



Q- l = 



Or 



1 

9 

2m 



dr / rfrlAI 



trln<T 



h -A 



(11) 



up to second order in fluctuations of /x(r, r) = fi + 
<fyt(r,r), h(r,r) = h + Sh(r,T) and A(r,r) = A + 
SA L (r, t) + i8A t (y 1 t) around their mean-field values. 
Here, we have introduced fi = (jx^ + /Uj,)/2 and h = 

mt - mi- 

By completing the squares in SA L and SA T , one can 
easily obtain: 









Tlhh) 








n^AT 



2? 



n 



A 1 /! 

A' 



n 



A 1 "?! 



(12) 



where we have introduced the polarization operators 
U ab (q,iu h ) = ^Yj 



tr r tr q+k)i6n+&l , h T)!,lT qjieti 



with r„ 



""3) T /i = 1) t a l — an d t a t = cr 2, and 



where the order parameter propagator is 



v- 1 = — 1 

9 



H^L A L IIa^ A r 

H A t A l IIa t A T 



(13) 



We omit the explicit expressions for these quantities, as 
they are straightforwar d g eneralisations of the expres- 
sions found e.g. in Rcf. [26| to the case h ^ 0. It is clear 
that Eq. (fl2|) is the generalisation of Eq. JU to the full 
matrix response and to nonzero q and e. 

The mean-field approximation to the mode spectrum 
is the solution of Eq. (fT0|) , with the response matrix given 
by the expression (fT2|) . Here, /j, and h are chosen so that 
n and m take the desired values, and inside the spinodal 
region, A is taken at the Sarma value corresponding to 
the maximum of fh(A, h) — as discussed, the stationary 
points of fh{A, h) coincide with those of f m (A, m). Evi- 
dently, a sufficient condition for a solution of Eq. (fl"Q|) is 
the occurrence of a pole at e(q) in the order parameter 
propagator, so we must solve det £> _1 (q, e(q)) = 0. In 
general, a numerical solution is called for. At the spin- 
odal lines, however, a diverging susceptibility implies a 
vanishing sound velocity, which we can find by expanding 
in e and q: 



A + P(e/q) iBe 
-iBe Qq 2 - Re 2 



(14) 
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l/k F a=0.21 

r\ Q<o 
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(b) 
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l/k p a=0.29 

Q>0 


(c) ' l/k F a=0.30 
Q>0 


(d) , 


l/k p a=0.77 
Q>0 
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FIG. 3: (Colour online) Upper panels: Unstable modes fre- 
quency uj/ef (w(q) = — ie(q)) versus momentum q/qF for dif- 
ferent values of the interaction strength 1/kpa across the spin- 
odal region and for fixed polarisation m/n = 0.5. Note that 
gaps correspond to complex frequencies. Lower panel: Plot 
of the most unstable mode frequency lo/ef (plus green), mo- 
mentum q/qF (times blue) and of the sound velocity 2c s /vf 
(solid black) across the spinodal region for fixed polarisation 
m/n = 0.5. 



where A = i£ k 6(£ k - h)£, B = i£ k e(£ k 



V SEl anCl 



8mA 2 



1 - 



3/2 3/2 



2Vh 2 



3/2 



A 2 

P(e/q) = W+L{e/v + q) + v_L{e/v.q)] 



with Ek = (e k — n) 2 + A 2 , and where v± = ^£^v± 
a,ndv± = dEk/dk\ e± are the density of states and velocity 
at the two solutions e± of E^ = h (take e± = = v± 
if there is no solution). Note that the phase stiffness 
Am A 2 Q is the supcrfluid density^, which changes sign 
inside the spinodal region (see Fig. |TJ) . 
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The higher order terms in the (1,1) entry of Eq. (fT4]) 
do not affect the sound velocity c s , which is the solution 
of 

[A + P{c s )] (Q - Rc 2 s ) - B 2 c 2 s = . (15) 

For obvious reasons, this analysis closely parallels that of 
Ref. for Bose- Fermi mixtures, and Eq. [15] reproduces 
the form of Eq. [5] found earlier. 

III. DISCUSSION 

At the spinodal lines A + P(0) = 0. This is the same 
condition that was used in Rcfs. Il6ll28l to identify the 
unstable region, although the possibility of metastabil- 
ity was ignored. A + P(0) < inside the spinodal re- 
gion, and one can distinguish two cases: When Q > 
the sound velocity c s is pure imaginary, while for Q < 
it is in general complex. As a consequence, the unsta- 
ble mode spectrum changes in character as one moves 
across Q = (see Fig. [3]). Approaching the region 
of negative superfluid density Q < from the region 
Q > 0, a second imaginary modes appears (see panel 
(c) in Fig. [3]) and the two merge (panel (b)): The re- 
sulting 'gap' region corresponds to complex frequencies, 
implying a 'flickering' component to the instability. The 
most unstable mode always corresponds to a pure imagi- 



nary frequency, however. While the most unstable mode 
frequency and wavevector go to zero on the superfluid 
side of the spinodal region, the instability towards a 
Fulde-Ferrell-Larkin-Ovchinnikov (FFLO) phase means 
that the characteristic wavevector does not go to zero on 
the normal side^. Indeed, one may view the early stages 
of spinodal decomposition as a transient FFLO state. 

The characteristic length and time scales at which in- 
homogeneities appear as the precursor of phase separa- 
tion are determined by the most unstable modes. At 
unitarity 1/kpa = and for Tp = lyuK, this time scale 
is of order of 400^s, which is on the same scale as the 
condensate formation time, as measured in Ref. [29l . The 
length scale is roughly 1/kp ~ 0.1/im. Both scales be- 
come larger as one approaches the spinodal lines, which 
will always occur somewhere in the trap. The unitarity 
region may be a suitable place to observe the late stages 
of spinodal decomposition and the possible existence of 
a coarsening regime. 

In conclusion, we have studied the early stage dynam- 
ics of phase separation in polarised Fermi supcrfluids. 
We expect that the investigation of these instabilities is 
within reach of current experiments. 
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One can in general show that, at the mean-field level, the 
free energy fT} is a function of m/n only, and conversely 
the grand canonical potential ([3]) is function of h/\fi\ only. 



